#ifndef ESTIMATOR_POLYA_CC
#define ESTIMATOR_POLYA_CC

#include<string>
#include<fstream>
#include<sstream>
#include<iostream>
#include<algorithm>
#include<iterator>
#include<vector>
#include<map>
#include<cstdlib>
#include<ctime>
#include<cmath>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_randist.h>
#include "def.hpp"
#include "est_polya_moment.cpp"
#include "est_polya_sgd.cpp"
#include "est_polya_newton.cpp"
#include "est_polya_fixedpoint.cpp"

using namespace std;

int EST_Polya(){
  // assume EST_BUFFER_INT is the input
  // we firstly perform a moment matchting method to get initial values
  // prepare random number generator
  gsl_rng *ge_gBaseRand;
  unsigned long ge_randSeed;

  ge_gBaseRand = gsl_rng_alloc(gsl_rng_mt19937);   /* specifying to use Mersenne twister MT-19937 as the uniform PRNG */
  ge_randSeed = rand();                    /* returns a non-negative integer */
  gsl_rng_set (ge_gBaseRand, ge_randSeed);    /* seed the PRNG */
  
  ofstream outFile;
  outFile.open(PuzaDEF::Instance()->output_file_name.c_str());

  double* estimated_alpha = new double[PuzaDEF::Instance()->EST_Polya_K];
  int* temp_n = new int[PuzaDEF::Instance()->EST_Samples];
  // clean
  for(int k=0; k < PuzaDEF::Instance()->EST_Polya_K; k++){
    estimated_alpha[k] = 0.0;
  }

  for(int i=0; i < PuzaDEF::Instance()->EST_Samples; i++){
    temp_n[i] = 0;
    for(int k=0; k < PuzaDEF::Instance()->EST_Polya_K; k++){
      temp_n[i] = temp_n[i] + PuzaDEF::Instance()->EST_BUFFER_INT[i][k];
    }
  }

  moment_matching(PuzaDEF::Instance()->EST_BUFFER_INT, temp_n,estimated_alpha, PuzaDEF::Instance()->EST_Polya_K, PuzaDEF::Instance()->EST_Samples);

  outFile << "[Moment Estimator]:" << endl;
  for(int k=0; k < PuzaDEF::Instance()->EST_Polya_K; k++){
    outFile << k << ":" << estimated_alpha[k] << " ";
  }
  outFile << endl;

  if(PuzaDEF::Instance()->currentPolyaEstimatorType == EST_POLYA_FIXEDPOINT){
    est_polya_fixedpoint(PuzaDEF::Instance()->EST_BUFFER_INT, temp_n, estimated_alpha, PuzaDEF::Instance()->EST_Polya_K, PuzaDEF::Instance()->EST_Samples);
    outFile << "[Fixed Point Estimator]:" << endl;
  }
  if(PuzaDEF::Instance()->currentPolyaEstimatorType == EST_POLYA_NEWTON){
    est_polya_newton(PuzaDEF::Instance()->EST_BUFFER_INT,temp_n, estimated_alpha, PuzaDEF::Instance()->EST_Polya_K, PuzaDEF::Instance()->EST_Samples);
    outFile << "[Newton Estimator]:" << endl;
  }
  if(PuzaDEF::Instance()->currentPolyaEstimatorType == EST_POLYA_SGD){
    est_polya_sgd(ge_gBaseRand,PuzaDEF::Instance()->EST_BUFFER_INT, temp_n, estimated_alpha, PuzaDEF::Instance()->EST_Polya_K, PuzaDEF::Instance()->EST_Samples);
    outFile << "[Stochastic Gradient Descent Estimator]:" << endl;
  }

  for(int k=0; k < PuzaDEF::Instance()->EST_Polya_K; k++){
    outFile << k << ":" << estimated_alpha[k] << " ";
  }
  outFile << endl;
  outFile.close();


  delete[] estimated_alpha;
  delete[] temp_n;  
  gsl_rng_free(ge_gBaseRand);
  return 0;
}

#endif
